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A general method is presented for estimating the uncertainty in hybrid models of gravitational waveforms 
from binary black-hole systems with arbitrary physical parameters, and thence the highest allowable initial 
orbital frequency for a numerical-relativity simulation such that the combined analytical and numerical 
waveform meets some minimum desired accuracy. The key strength of this estimate is that no prior 
numerical simulation in the relevant region of parameter space is needed, which means that these techniques 
can be used to direct future work. The method is demonstrated for a selection of extreme physical 
parameters. It is shown that optimal initial orbital frequencies depend roughly linearly on the mass of 
the binary, and therefore useful accuracy criteria must depend explicitly on the mass. The results indicate 
that accurate estimation of the parameters of stellar-mass black-hole binaries in Advanced LIGO data or 
calibration of waveforms for detection will require much longer numerical simulations than are currently 
available or more accurate post-Newtonian approximations — or both — especially for comparable-mass 
systems with high spin. 



Numerical relativists face a thorny dilemma when creat- 
ing initial data for simulations of binary black holes. Two 
competing motivations vie for control of one key parameter: 
the initial orbital frequency of the binary, Q - On one 
hand, the larger the initial frequency is, the more quickly 
the simulation will run. In fact, at lowest order, doubling 
Qq will shorten a simulation by a factor of six. On the other 
hand, the smaller the initial frequency is, the clearer the 
correspondence will be between the numerical simulation 
and the system found in nature. Post-Newtonian approxima- 
tions will be more accurate; the velocity and spin of each 
black hole will be more well defined and easily measured; 
even the junk radiation will be smaller [1]. In practice, 
numerical relativists choose Qq largely by intuition, with 
primary considerations being the available computer time 
and the roundness of the number. This lack of precision can 
lead to simulations that are too short and should ideally be 
redone, or are longer than necessary — a waste of resources 
in either case. More objective choices are possible, and 
will be needed to improve the effectiveness of numerical 
relativity. This paper demonstrates a technique 1 to estimate 
the optimal value of Q , even for systems in unexplored 
regions of parameter space. 

The field of numerical relativity (NR) exists because of 
the failure of post-Newtonian models (PN); at some point 
NR must take over from PN. Of course, PN approximations 
don't simply break down at one catastrophic instant, having 
been perfectly accurate before. Rather, the approximations 



1 The technique expands on one introduced in Ref. [2] — which was 
partially implemented in Ref. [3] — to apply to the time-domain methods 
currently in use by most numerical relativists, to use the most accu- 
rate inspiral models available, and to include more general accuracy 
requirements. The results obtained here broadly agree with the results 
of Refs. [4-6], which test whether completed numerical simulations are 
long enough. 



gradually deteriorate as they approach merger. The ques- 
tion of exactly where NR needs to replace PN is thus a 
question of how accurate the model needs to be. In the 
context of designing model waveforms for detection, we are 
given precise objectives. This allows us to quantitatively 
resolve the conflict between decreasing the length of a 
simulation and improving the quality of the final modeled 
waveform. The quality of the final waveform is impacted 
by the accuracy of both PN and NR data. Because we 
already have the PN data that will be used in the final 
waveform, we can test how much of it can be used if we 
are to achieve a target accuracy. Where that ends is where 
the NR simulation must begin. 

Estimating the impact of PN errors depends on under- 
standing how the data will be used in the finished product. 
In this context, that means understanding waveforms used 
in data analysis for gravitational-wave detectors. Advanced 
detectors of the near future will be more sensitive over 
broader ranges of frequencies than current detectors, which 
tightens the requirements for accurate modeling of physical 
waveforms [7, 8]. In particular, analysis of detector data 
will require waveforms that are not only more precisely 
coherent, but also coherent over a greater range of frequen- 
cies. PN waveforms alone are expected to be sufficient for 
detection up to a total system mass somewhere 2 around 
12 M Q [9], and are expected to be entirely inadequate for 
parameter estimation. Exclusively numerical waveforms, 
on the other hand, are limited in their usefulness to high- 
mass systems; current numerical simulations can only cover 
the Advanced LIGO band down to 10 Hz for masses greater 



2 The exact mass delimiting the range of validity depends, of course, on 
the parameters of the particular system and the detector in question. 
Moreover, this says nothing about parameter estimation. The methods 
of this paper are primarily applicable to parameter estimation, but do 
have some relevance to detection. 
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than about 100M G [10]. This leaves a large gap in an 
astrophysically interesting range [11]. To close the gap 
using NR alone would require simulations roughly 250 
times longer than are currently available — or more for 
parameter estimation. Notwithstanding possible dramatic 
improvements to the efficiency of long simulations [12-14], 
this will be impractical for some time to come, especially 
for large surveys of parameter space [15]. 

Instead of running such dramatically long simulations, 
we "hybridize" — synthesizing a single coherent waveform 
by combining the long PN inspiral with the short NR merger 
and ringdown. In the interest of simplicity of presentation, 
let us assume for the moment that a hybrid waveform is con- 
structed by aligning the PN and NR waveforms at a single 
point in time — using PN data before that point and NR data 
after — and that the NR portion is essentially perfect. Now, 
as the alignment point moves closer to merger, the hybrid 
waveform will become less accurate since it includes more 
of the deteriorating PN waveform, as demonstrated in Fig. 1 . 
In fact, because the hybridization procedure uses only the 
worst data when aligning at very late times, errors in the PN 
waveform just prior to the alignment point will unavoidably 
taint the data, meaning that the long inspiral and the merger 
will be dephased relative to the physically correct waveform. 
The error in the complete hybrid, then, depends crucially on 
the error in the PN waveform and particularly on the growth 
in the error at late times. 

We cannot know the error in a PN waveform — being 
the difference from the unknown correct waveform it at- 
tempts to model — without unreasonably long numerical 
simulations. We can, however, estimate our uncertainty 
by creating a range of equally plausible PN waveforms. 
We can then attach these to some ersatz NR waveform, 3 
such as effective one-body (EOB) or phenomenological 
waveforms [16-22], to form a range of plausible hybrids. 
We can compare that range quantitatively to the given error 
budget. By repeating this process using each of many 
possible hybridization frequencies, we can discover which 
hybridization frequencies produce waveforms that satisfy 
the error budget. The highest such frequency minimizes 
the length of the simulation, and is thus the optimal value. 

To summarize, the work to be done in this scheme 
consists of three basic steps: 

1. Choose an appropriate ersatz NR waveform and con- 
struct a range of plausible inspiral waveforms (Sec. I); 

2. Hybridize the inspiral and ersatz NR waveforms, evalu- 
ate the mismatches among them, and repeat for various 
hybridization frequencies (Sec. II); 

3. Choose Q.q to agree with the highest hybridization fre- 
quency that satisfies the error requirements (Sec. III). 
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3 As shown later in this paper, the final results will not depend strongly 
on the particular choice of ersatz NR waveform. 
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FIG. 1. Error caused by aligning at late times. This plot 
shows the phase error in hybrid waveforms created by aligning 
an ersatz NR waveform (EOB) to a PN model at various times, 
relative to the hybrid when aligned very far in the past. As the 
alignment point is moved closer to merger (t/M = 0), the total 
phase error increases because the hybrid waveform incorporates 
more and more of the inaccurate PN data. Note that, in each case, 
the error grows most rapidly near merger. If the merger occurs at 
a frequency to which the detector is sensitive, this phase error will 
negatively impact the match found by data analysts (see Sec. I A). 
Evaluating that impact using artificial NR data is the essence of 
the method presented in this paper. 



This, of course, assumes that there is some freedom in 
choosing the initial orbital frequency, as will be the case for 
NR groups setting out to run a simulation; there is always a 
choice between running many short simulations and fewer 
long simulations — or requesting a larger allocation of com- 
puter time. This method is designed to help in making that 
choice. On the other hand, given a completed waveform, 
steps 1 and 2 above can be used to evaluate the uncertainty 
in the resulting hybrids. Finally, understanding the results 
can aid in the design of reasonable and effective accuracy 
goals, before any simulation is undertaken. 

The method will be demonstrated for a few interesting 
cases, probing the "corners" of a simple parameter space: 
equal-mass nonspinning, equal-mass high-spin, large mass- 
ratio nonspinning, and large mass-ratio high-spin systems. 
The key result of this paper is the plot of the uncertainty of 
the hybrids for those systems, Fig. 4, discussed in Sec. Ill C. 
Finally, Sec. IV summarizes the conclusions and outlines 
possible applications and extensions to this method. 

The general method presented here can be applied quite 
broadly. However, to demonstrate the method, we need to 
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make several specific choices. Details of the implementa- 
tion are given in Appendix B. In particular, this paper uses 
a simplified EOB model, described in Appendix B 2, to 
supply the ersatz NR waveforms. In Appendix A, the results 
for the equal-mass nonspinning case are redone using actual 
numerical data, to check — to the extent possible — that the 
final results are not sensitive to this choice for the ersatz 
NR waveform. The uncertainty will be measured by a 
range of plausible waveforms formed by hybridizing the 
EOB merger and ringdown with the inspiral portion of EOB 
and TaylorTl-T4 approximants [23-25] using all known 
information (full PN orders), as recently recalculated and 
summarized by members of the NINJA-2 collaboration [26]. 
Specifically, the amplitude includes terms up to 3.0 PN 
order, and the phase includes terms up to 3.5 PN order. See 
Appendix B 1 for more detail. The hybridization will be 
done in the time domain by aligning at particular frequen- 
cies [25, 27], then blending the PN and NR waveforms 
as described in Ref. [10]. The resulting hybrids will be 
compared along the positive z axis, computing the match 
[see Eq. (3)] using the Advanced LIGO high-power noise 
curve with no detuning [8], scaling the total system mass 
between 5 M Q and 50 M Q . 

Throughout this paper, the uppercase Greek letters <D 
and Q refer to the orbit of a binary, contrasting with the 
lowercase Greek letters (f> and u> which refer to the phase 
and frequency of the emitted gravitational waves. Unless 
otherwise specified, <f> and u> refer to the (£, m) = (2, 2) 
mode in a spin s = — 2 spherical harmonic decomposition 
of the gravitational wave. 

I. CREATING THE WAVEFORMS 

The first task before us is to construct a large group 
of model waveforms to be compared to each other. In 
this section, a review of the standard error measure used 
in data analysis motivates the use of complete inspiral- 
merger-ringdown waveforms for our tests, while finding 
encouraging signs that high accuracy of these waveforms 
is not essential. We then examine in greater detail the 
construction of a credibly broad selection of PN waveforms. 

A. Motivation 

In constructing of a range of plausible model waveforms, 
we need to understand the ultimate form of measurement 
when designing templates for gravitational-wave detection: 
the match. This quantity is based on the inner product 
between two waveforms defined as the integral of the 
noise-weighted product of the signals in the frequency 
domain [28]: 
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where S„(f) is the power spectral density (PSD) of noise in 
the detector, and 8<p(f) is the phase difference between the 
two waveforms in the frequency domain. In the following, 
we will assume / > 0. The second form shown here demon- 
strates a useful way to understand the inner product, by 
separating the integrand into two factors. We begin with the 
first factor, the ratio of amplitudes to noise. The amplitude 
for an example system and the Advanced LIGO noise curve 
are plotted in Fig. 2. The height of the amplitude above 
the noise curve shows when the inner product can rapidly 
accumulate; more height means more rapid contribution to 
the inner product. Of course, rapid contributions cannot 
result in a large inner product unless those contributions 
are also coherent. This requires the second factor of the 
integrand in Eq. (lb) to be constant. In particular, a large 
inner product requires dip to be very close to zero across the 
full range of frequencies for which the signal amplitude is 
significantly larger than the noise. 

While the amplitude of the waveform from a given sys- 
tem is fixed, 6(f> has two inherent degrees of freedom we 
can adjust to improve the inner product. These degrees of 
freedom are related to the fact that the merger time and 
orientation of an astrophysical binary are unknowns that 
must simply be measured. Specifically, we are free to shift 
the time and phase of either time-domain waveform by AT 
and AO. In the frequency domain, this has no effect on the 
amplitude of the waveform (the curves of Fig. 2 will not be 
affected), but 6(f> changes roughly as 



6<Kf)->6<Kf) + 2A<l> + 27rfAT, 



(2) 



We can use this time- and phase- shift freedom to ensure 
that the phase difference between two waveforms is smallest 
at frequencies for which the detector is most sensitive to 
that particular system, and thus maximize the inner product. 
The maximum possible (normalized) inner product [28] is 
called the match, which will be our basic measure of error: 



(h a \h b ) = max 



(h a I h h ) 



AT, A* yj(h a \h a )(h b \h b ) 



(3) 



This quantity takes a value between (for completely 
dissimilar waveforms) and 1 (for identical waveforms). Be- 
cause many of the matches encountered below will be very 
close to 1 , it is preferable to use another quantity called the 
mismatch [29-31], which is given by 



MM(h a ,h b ) := l-(h a \h b ) 



(4) 



Here, values close to indicate the waveforms are simi- 
lar. The maximum possible signal-to-noise ratio (SNR) at 
which a given signal s can be detected is given by 



(5) 



The mismatch MM (s, h) between a signal s and a template 
h is essentially the percentage loss in SNR due to errors in 
the template [32]. 
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FIG. 2. Waveforms in the frequency domain. Amplitudes 
from an equal-mass nonspinning binary are shown, scaled to total 
system masses of 10 M @ and 100 M Q at 100 Mpc, and compared 
to the noise spectral density of Advanced LIGO (both quantities 
in units of strain-Hz -1 ^ 2 ). The factors multiplying the amplitudes 
are chosen to account for the logarithmic scaling of the horizontal 
axis, the factor of 2 in Eq. (1), and the fact that only positive 
frequencies are plotted. The triangles on the waveforms show 
the approximate initial frequency of the longest numerical simu- 
lation currently available. The circles show the frequency of the 
innermost stable circular orbit (ISCO) — basically the frequency 
at which PN approximations are expected to be useless. This plot 
shows that, for the 10M Q system, nearly half the contribution to 
the inner product of Eq. (1) comes from the PN data (to the left 
of the triangle), and the rest from the NR data (to the right). For 
the 100 M e system, on the other hand, the inner product is given 
almost exclusively by NR data. 



This understanding of the match teaches us two very im- 
portant lessons that will guide our approach to our problem. 
First, the optimal phase shift — and thus the value of the 
match — will depend on the relative distribution of power 
at different frequency bands. If, for example, our model 
waveform simply ends at ISCO, it will fail to model the 
a substantial portion of the physical waveform. In that 
case, the maximization of Eq. (3) will not need to balance 
the dephasing between the two portions of the waveform, 
for example. The requirements for phase coherence will 
be much looser than they, in fact, need to be. Thus, our 
model waveforms must have roughly the same distribution 
of amplitude as the physical waveform, across the entire 
sensitive frequency band of the detector. If our objective 
is to evaluate matches before any numerical simulation 
is done, we will need a suitable approximation to the 
merger/ringdown waveform. For that purpose, this paper 



uses the EOB waveform [16-19], which extends through 
merger to ringdown. Other complete waveforms could also 
be used [20-22]. Fortunately, we will see that the final 
results will not depend strongly on the choice of ersatz NR 
waveform. For example, after a simulation is done, we can 
go back and check that the results agree if we use the NR 
waveform itself. This is done in Appendix A for the equal- 
mass nonspinning system. Even more extreme, the orig- 
inal proof-of-principle demonstration of this method [2] 
used stationary-phase approximated (SPA) waveforms ter- 
minated at the light-ring frequency. For small mismatches, 
the results achieved in that test compare well to the results 
achieved in this paper, even though the SPA waveform 
mismatches the numerical and EOB waveforms by more 
than 8 % over the relevant mass range [10]. 

The second lesson gleaned from these considerations of 
the match is that the phase of a model waveform does not 
come into the match; only the phase difference between 
models matters, as shown explicitly in Eq. (lb). This fine 
distinction has real importance for us because it implies 
that the phase error in our ersatz NR waveform relative to 
the correct physical waveform is not important. Certainly 
the final model waveform should resemble the physical 
waveform as closely as possible, but we will assume that 
errors in the portion of the final waveform covered by NR 
data will be accounted for separately in the error budget — 
or are essentially negligible compared to the PN errors. 4 
When comparing two plausible waveforms hybridized with 
the same ersatz NR data at some frequency /hyb, the phase 
difference during the NR portion (f > / n yb) will De zero to 
a very good approximation, at least until the waveforms are 
shifted in time and phase to maximize the match. But even 
then, the phase difference will not depend in any way on the 
phase of the ersatz NR waveform; by Eq. (2), it will be 



<W) = 2A<D + 27r/Ar for/>/ t 



hyb- 



(6) 



Thus, the mismatch between plausible hybrids is not di- 
rectly sensitive to the particular phasing of the ersatz NR 
waveform. Of course, that phasing will affect the alignment 
during hybridization, which can affect the relative power in 
different portions of the waveform or the function 6(/>(f) for 
frequencies / < /hyb- However, the results below show 
that the ersatz waveform is not dominating the uncertainty, 
which suggests that even this phase error is not important. 

These considerations lead us to conclude that the ersatz 
NR waveform must be reasonably accurate, especially in 
terms of modeling the relative power in different parts of 
the waveform. However, we are also given hope that our 



4 Of course, this method readily applies to quite general error budgets. 
For example, if the expected uncertainty in the numerical waveform can 
be estimated, even if that error depends on the initial frequency of the 
simulation, the error budget can be trivially extended to include that 
estimate. See Sec. Ill E for more details. 
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final result — the predicted value for — does not depend 
strongly on the accuracy of that ersatz NR data. We can test 
this expectation and will see in Appendix A that, at least in 
the case of the equal-mass nonspinning system, it is indeed 
well founded. 



B. The range of plausible PN waveforms 

Reliable final results depend on an accurate range of 
plausible hybrids, correctly depicting the possible error in 
the analytical models. The range must be neither too broad 
nor too narrow: too broad, and we will conclude that the 
error is large, and thus begin the simulation earlier than 
necessary; too narrow, and we will be overconfident in our 
models, and thus waste time producing an inaccurate hybrid 
waveform. Unfortunately, we have no way of knowing the 
error in our analytical models before the fact. If, however, 
we assume that our models are not wrong, but are simply 
incomplete, we should be able to trust the uncertainties in 
the model to estimate the error. 

Analytical relativity has produced multiple methods of 
calculating waveforms for black-hole binaries. Roughly 
speaking, these different methods should be equivalent at 
the level of our knowledge of the true waveform. To 
the extent that they are different, they are uncertain. In 
fact, we will use precisely this range of differences as our 
range of uncertainty, and thus our estimate for the errors 
in the analytical waveforms. Thus, choosing our range of 
plausible hybrids comes down to choosing representatives 
of the various methods for calculating analytical waveforms. 
The representative methods of calculation we will use are 
the TaylorTl-T4 [23-25] andEOB [16, 33] models. 

An objection might be raised that the EOB waveform 
is more accurate and, in particular, "breaks down" more 
elegantly than the basic PN waveforms; the implication 
being that the TaylorTn waveforms should not be included. 
It would certainly be possible to use EOB waveforms alone, 
employing so-called flexibility parameters [34] to delimit a 
range of plausible waveforms, for example. Unfortunately, 
while EOB waveforms can be tuned very precisely to re- 
semble the late-time behavior of numerical waveforms after 
the fact [33, 35], there is no evidence that the inspiral — 
which is of more interest here — will be more accurate 
after this tuning [36], or that any portion of the EOB 
waveforms will be more accurate before such tuning can 
be done. In fact, Blanchet has suggested that EOB appears 
to converge toward a theory which is different from general 
relativity [37]. Because the results will need to apply in 
regions of parameter space where no simulation has yet 
been done, and because this paper attempts to reflect meth- 
ods currently in use by the community [26], we will take 
the more conservative approach of including the TaylorTn 
approximants. 



II. EVALUATING UNCERTAINTY 

With a selection of plausible waveforms in hand, we can 
now evaluate the differences between them in terms of the 
match. The aim is to treat the artificial data (the ersatz 
NR waveform) exactly as data from the actual simulation 
will be treated. First, methods of hybridization will be 
reviewed. This hybridization will be performed for a variety 
of hybridization frequencies, <Dh y b- Then, we can simply 
evaluate the match between the various hybrids at each tUh y b, 
as a function of the total mass of the system. 

A. Hybridization techniques 

Combining inspiral and merger/ringdown waveforms is a 
delicate process, beginning with the procedure for aligning 
the waveforms by matching up the arbitrary time and phase 
offsets in the data. As described by MacDonald et al. [6], 
this part of the process has large potential effects on the 
accuracy of the final result; in their example a misalignment 
of just 1 M in the time values of the two waveforms resulted 
in a mismatch of up to MM = 0.01. Clearly, this part of the 
process must be handled carefully. Many techniques have 
been devised for doing so, resulting in a variety of choices 
to be made. 

First, alignment of the time and phase offsets may be done 
in either the time domain or the frequency domain. For 
the particular case of hybridizing to numerical waveforms, 
the numerical data is short, beginning at high frequency 
and — in particular — having large amplitude. Transforming 
such data into the time domain will either introduce Gibbs 
phenomena, which will spoil much of the NR data, or 
require windowing, which will waste much of the NR data. 
Therefore, we use time-domain alignment for our purposes, 
as is used throughout most of the current literature/ 

Second, we must choose a criterion for deciding how 
well the two waveforms are aligned after offsetting the 
time and phase. Many possibilities have been suggested 
for this purpose, including the magnitude of the differ- 
ence in the complex h data [20]; the gravitational- wave 
phase and frequency [25]; even the orbital phase and fre- 
quency [24]. The gravitational-wave phase and frequency 
are often chosen by numerical-relativity groups to produce 
hybrids because the post-Newtonian phase is known to 
higher order than the amplitude, and is thus more likely to 
result in an accurate alignment. Because of its popularity 
and simplicity, phase alignment is used in this paper. 

Finally, the alignment procedure depends on the width 
of the region over which the criterion chosen above is 
evaluated. For example, to align the phase (and implicitly 



5 In the original proof-of-principle demonstration of the method de- 
scribed in this paper [2], frequency-domain alignment was used. In 
that case, the ersatz NR and PN waveforms were various versions of 
the TaylorF2 waveform — which is calculated in the frequency domain, 
which means Gibbs phenomena are not relevant. See also Ref. [4]. 
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frequency) of two waveforms, a common method [10] is to 
minimize the squared difference between them: 



E(AT, A0) 



- f [0: 
Jt, 



nr(0 - 4>PN(t + AT) - A<p] 2 dt . 



(7) 

The alignment then depends on t\ and ti- There is a certain 
trade-off here, between using a short region so that less of 
the inaccurate PN waveform can be used, and using a long 
region to smooth out any irregularities in the numerical 
data, such as junk radiation or residual eccentricity. Ad- 
ditionally, the range \t\, tj\ must capture some curvature in 
the graph of §{t) for an accurate alignment, which means 
that the range must become larger at lower frequencies. Ref- 
erence [6] suggests a simple but robust method of choosing 
this range, where t\ and t% extend to frequencies 5 % above 
and below some central frequency. This ensures that the 
range is neither too large at high frequencies, nor too small 
at low frequencies. 

In our case, the "numerical" data is the EOB waveform, 
which has essentially no eccentricity or noise. Thus, for 
simplicity of presentation and implementation, we take the 
limit of this procedure as ti approaches t\ . 6 To do so stably, 
we adjust the time offset AT so that the frequencies are the 
same at some time t^, then adjust the phase offset Ac/) so 
that the phases are the same at fhyt> : 

^NR<Ayb) - WpN^hyb + AT) , (8a) 
0NR(?hyb) = ^PN^hyb + AT) + A<p . (8b) 

Here, the optimal offsets will depend on the time at which 
the alignment condition is imposed. Below, this depen- 
dence is described using the frequency itself: Whyb := 

<^NR(?hyb)- 

The alignment just described is typically applied to the 
(t, m) = (2, 2) mode in a spin-weighted spherical harmonic 
decomposition of the gravitational waves. Other modes 
must also be aligned. However, we have fixed the only 
degrees of freedom, which means that the other modes are 
already determined. In general, the amplitude and phase of 
any mode of the PN waveform is transformed according to 



A ( M -» A%(t - AT) , 



(9a) 

0pn W -» 0pn (' - AT) - m A0/2 + 2 n n , (9b) 



for some integer n that ensures continuity of the phase. 

Now, having aligned the two waveforms, we need to 
produce a single waveform. Because only <p 2 ' 2 is guaran- 
teed to be continuous, discontinuities are possible in other 



6 It has been checked that the results of this paper are essentially identical 
when using Eq. (7) with values of t\ and tj as prescribed in Ref. [6]. 
This is a better choice for data from simulations, and is more typical 
of hybridization as practiced by numerical-relativity groups, but would 
introduce an unnecessary layer of complexity to the discussion here. 



quantities, so a hybrid is usually formed using a transition 
function to blend the two waveforms. Here, because the dis- 
continuities are mild, we use a basic linear transition [ 10] of 
width 10 M centered on the alignment point, for amplitudes 
and phases of all modes. 

A special case arises when the PN waveform is the EOB 
approximant. Rather than actually splitting the EOB wave- 
form into two parts and recombining them, the complete 
EOB waveform is used as the EOB "hybrid". 

B. Mismatches 

Now, having formed a series of hybrids using various 
PN approximants, and a range of hybridization frequencies 
Whyb: we can evaluate the difference between them using 
the same criterion as will be used with the numerical data. 
In particular, we take the mismatch [Eq. (4)] using the 
Advanced LIGO zero-detuning, high-power noise curve [8]. 
The waveforms are projected onto the positive z axis us- 
ing all available modes of the spin-weighted spherical- 
harmonic decomposition. Post-Newtonian calculations 
have been carried out through t = 8: 



h(t) = % 



(=2 m=-l 



A e '"\t)^ 



- 2 r, m (o,o) 



. (10) 



Note, however, that the quasinormal-mode portion of 
the EOB waveform has only been extended to include 
(£,m) e {(2, ±2), (2, ±1), (3, ±3), (3, ±2), (4, ±4)} [33]. 
Other modes are set to zero during ringdown. 

As a first example, Fig. 3 shows the mismatches between 
each pair of hybrids for the equal-mass nonspinning system 
scaled to a total mass of 20 M Q . At any particular value 
of o>h y b there is a range of mismatches, indicating that 
some pairs of hybrids happen to agree with each other very 
closely, while some are quite different. There is no reason to 
suspect that a pair of hybrids in close agreement with each 
other also agree with the exact waveform. Rather, if these 
are all plausible waveforms, the uncertainty in our model is 
given by the maximum mismatch between any pair. For this 
particular system, that pair happens to be the TaylorTl and 
T3 waveforms for most values of o»h y b, though in general 
TaylorTl is most dissimilar from the other waveforms. 

We can follow the maximum mismatch as a function of 
frequency, and notice a general trend: increasing hybridiza- 
tion frequency results in increasing uncertainty. This is to 
be expected for two reasons. First, the PN approximation 
should be very accurate at low frequencies, but break down 
at higher frequencies. For example, the time and phase at 
which o>h y b occurs in the PN waveform will become more 
uncertain as o»h y b increases. This results in uncertainty in 
the alignment between the two parts of the hybrid. 

Second, as the hybridization frequency increases, the 
detector will simply be more sensitive to the differences 
between hybrids. The upper horizontal axis of Fig. 3 shows 
the physical hybridization frequency /hyb- Comparison 
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/hyb (Hz) 




FIG. 3. Mismatches between hybrids as a function of w hy b 
This plot shows the mismatch between pairs of hybrids using 
different approximants, for the equal-mass nonspinning system 
with M lot = 20 M Q . At any particular value of co^, the maximum 
mismatch between each pair of hybrids is the uncertainty in the 
final waveform hybridized at that frequency. If our target accuracy 
were, for example, MM < 10" 3 for this system mass, this plot 
shows us that the NR waveform would need to contain the GW 
frequency of G M lot Whyb / c 3 = 0.02, which naturally implies the 
initial orbital frequency of the simulation, Q.q. Achieving a smaller 
uncertainty requires lower <x>hyb, corresponding to times at which 
the PN approximants agree more closely. 



with the noise curve plotted in Fig. 2 demonstrates that the 
mismatch grows very quickly just as the hybridization point 
passes the "seismic wall" of the detector (/ se ismic — 10 Hz) 
where the sensitivity is improving rapidly with increasing 
frequency. The mismatch then begins to level out as the 
detector sensitivity levels out. 

In the next section, we discuss how to use plots like 
the one in Fig. 3, along with a target accuracy, to find the 
optimal initial orbital frequency of a numerical simulation. 

Before moving on, however, let us pause to note an 
interesting feature of the last plot. The various comparisons 
separate into distinct groups. The largest mismatches in- 
volve TaylorTl hybrids (solid lines), and all hybrids using 
Tl have large mismatches with other hybrids. Setting aside 
the Tl waveforms, we see a similar trend develop at high 
frequencies for T3 hybrids compared to hybrids other than 
Tl (dashed lines). We could even push this categorization 
to the T2 hybrids compared to waveforms other than Tl 
and T3 (dotted lines), though this leaves little room for 
comparison. We might say that Tl is dominating the un- 
certainty, in the sense that it is furthest from the consensus 



of other hybrids. At higher frequencies, T3 departs from 
that consensus, followed by T2, leaving only the T4 and 
EOB agreeing with each other. What is striking about this 
pattern of disagreements is that it is identical to the pattern 
in errors relative to the numerical waveform [25], where 
Tl is least accurate, followed by T3, with T2 slightly less 
accurate than T4; the T4 and EOB waveforms agree with 
the numerical result nearly within numerical errors. 

An optimist might suggest that close systematic agree- 
ment between two models like T4 and EOB is unlikely — 
given the size of the function space through which they 
are free to roam — unless they also agree with the exact 
waveform. This would imply that we should take the 
small mismatch between T4 and EOB as the uncertainty in 
those waveforms. Or, slightly less optimistically, we might 
discount Tl as being too far from the other waveforms, and 
thus a mere anomaly. Unfortunately, similar patterns do not 
develop for the other systems investigated below. For now, 
we leave this as a mere observation, and take the largest 
mismatch as an indicator of the uncertainty in any given 
model. 

III. USING MISMATCH TO FIND Q 

Given a target accuracy, the uncertainty implied by Fig. 3 
suggests a natural starting point for the numerical simu- 
lation of that system, simply because the simulation must 
include — at a minimum — the corresponding a)^. m tn i s 
section, the uncertainty estimate of the previous section 
is used to produce an optimal initial orbital frequency for 
the simulation. This is then generalized to apply across a 
range of masses and to incorporate more complicated target 
mismatches. In the process, the uncertainties for a small 
selection of astrophysical systems are shown. 

A. Optimizing Dq for a particular mass 

Figure 3 establishes a relationship between the uncer- 
tainty in plausible hybrids and the frequency at which they 
are hybridized, Whyb- If, f° r example, we wish to model an 
equal-mass nonspinning system of total mass M tot /M - 
20 with a target accuracy of MM target < 10~ 3 , this plot 
demonstrates that the final hybrid waveform must be formed 
with GMtot^hyb/c 3 ^ 0.02. Naturally, the numerical 
simulation must include that frequency. There are two 
simple methods for turning the GW frequency into an initial 
orbital frequency for the simulation, Qq. 

First, we might use the basic approximation Qq ~ Whyb/2. 
In this case, the result above would suggest an initial 
frequency of G M tot flo/c 3 < 0.01. 7 The actual simulation 



7 This is roughly half the initial frequency of current long simulations. To 
lowest order, the length of a simulation goes as T oc d g 8 ' 3 , which means 
that a simulation held to this standard needs to be roughly 2 8 ' 3 ss 6.3 
times longer than current long simulations. 
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should probably begin somewhat earlier than this, to allow 
junk radiation to leave the system, and to ensure that the 
alignment region of Eq. (7) does not extend past o»hyb- 
These considerations depend on the particular formulation 
of Einstein's equations and numerical methods used in the 
simulation, and are thus beyond our scope. Ultimately, the 
numerical relativist will use his or her judgment to produce 
some time At before Wh y b/2 at which the simulation should 
begin. In this regard, an additional PN approximation may 
be useful: 

co hyh /w h yb\ 11/3 96v /GM,„,f 



where v = ni\ mil(m\ +m2) 2 is the symmetric mass ratio of 
the individual black holes. The extra term is derived from 
the lowest-order PN approximation for the evolution of the 
orbital frequency [37]. Because of the approximations, this 
method may fail in certain extreme cases. 

Alternatively, and more robustly, we might simply refer 
to any of the PN models contributing to our estimate, to find 
the orbital frequency corresponding to a>hyb- Moreover, if 
some At is prescribed for the simulation, the PN model can 
be used to find the orbital frequency Qq occurring at a time 
At before the GW frequency a>hyb- 

The particular example just discussed applies only when 
the target accuracy is required for M tot - 20 M Q . This 
would be relevant to the situation where, for example, a 
source has been detected, and its parameters are known 
to reasonable accuracy, but further simulations are being 
done for accurate parameter estimation. More generally, 
however, we should expect to encounter broader accu- 
racy requirements, which might apply across a range of 
masses [15]. The rest of this section will extend this 
example to account for various masses; to demonstrate the 
uncertainty for a selection of interesting systems; then to 
allow the target accuracy to vary as a function of mass; and 
finally to allow the target accuracy to vary as a function of 
both mass and hybridization frequency. 

B. Optimizing Oo for a range of masses 

The mismatch curves plotted in Fig. 3 depend strongly 
on our choice of the total system mass. To generalize this 
to be a function of both Whyb an d M tot , we create the contour 
plot in the upper left of Fig. 4. Slicing through that plot at 
A^tot/-^o - 20 gives the uppermost curve in Fig. 3. For 
comparison, this quantity is also plotted for a selection of 
systems with different mass ratios or spins, as discussed in 
greater detail in Sec. IIIC. 

Again, given an accuracy requirement, we can use this 
plot to derive the optimal initial orbital frequency. If the 
requirement is a target accuracy of MM target < 10~ 3 , we can 
follow the 10~ 3 contour in the plot, and see that it is always 
above G M tot a>h y b/c 3 ~ 0.0075 for the range of masses 
shown. As before, the initial orbital frequency Oo is then 
deduced from this value by using Eq. (1 1) or by consulting 



the PN model, as described above. This stringent accuracy 
requirement calls for numerical simulations roughly 87 
times longer than the longest current simulation of this 
system. 

C. Comparing the uncertainty in various systems 

While the equal-mass nonspinning case nicely illustrates 
the method of finding optimal initial frequencies, systems at 
the boundaries of current numerical capabilities also hold 
a great deal of interest. Fig. 4 illustrates the uncertainty for 
four systems: 

1. Equal-mass, nonspinning; 

2. Equal-mass, aligned spins^fi = Xi = 0.95; 

3. Mass ratio 10:1, nonspinning; 

4. Mass ratio 10:1, aligned spins x\ — Xi — 0.95. 

The quantities x\ an d X2 are the components of the di- 
mensionless spins along the orbital angular- momentum 
vector. In the lower-left plot, case (3), frequencies above 
GM to tWhyb/c 3 ~ 0.055 are not possible because the 
TaylorT3 approximant ends at that frequency for this system. 
We can treat those higher frequencies as having MM = 1 . 
Note that the smallest known black holes have masses of at 
least M « 3 M Q [38]. This suggests that the smallest total 
mass for systems with q = 10 would be M tot ~ 33 M Q . 
Thus, the low-mass regions of those two plots may not be 
interesting astrophysically. 

In each case, we see the basic trends noted in Sec. II B. 
Two factors drive the mismatch: how far Whyb has entered 
the sensitive band of the detector, and how poorly the post- 
Newtonian approximation performs up to that frequency. 
Each plot of Fig. 4 includes a dotted red line denoting 
/seismic the lower bound of sensitivity in Advanced LIGO. 
Below this line, the mismatch must be zero, because the 
data in the detector's sensitive band is identical for any two 
waveforms — it is just the ersatz NR data. As we move above 
this line, a larger fraction of the data in the corresponding 
hybrids comes from different approximants. Thus, the 
mismatch increases. We expect it to increase more quickly, 
as a function of Whyb f° r systems that are not well described 
by PN approximations. Indeed, comparing the plots, we see 
that any given contour line moves closer to /seismic as either 
the mass ratio or spin parameter increases. This is simply 
an explicit confirmation that the physical system is not well 
modeled in more extreme cases. 

Though less important for our purposes, the contour 
lines for larger mismatches also obey a similar bound. In 
each plot, the dashed green line shows f2o%, the twentieth 
percentile of power in the match. That is, assuming 5<p = 
in Eq. (lb), f%o% is the frequency for which the integral has 
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FIG. 4. Maximum mismatch between plausible hybrids for a selection of systems. These plots show the maximum mismatch between 
any pair of hybrids formed by hybridizing TaylorTl-T4 or EOB with an ersatz NR waveform (EOB) at a^yh, as described in Sec. II A, 
when scaled to a total system mass of M tot . This quantity describes the uncertainty in the models used to form the hybrid. The plots show 
data from four cases with mass ratios denoted by q, which is the ratio of the larger to smaller mass, and components of the dimensionless 
spins x aligned with the orbital angular-momentum vector. Note that M lot < 33 M G may not be interesting astrophysically for black-hole 
binaries when q = 10. The region G M^oj^/c 3 > 0.055 is inaccessible in the q = 10, x = case, because the TaylorT3 approximant 
ends at that frequency; this could be interpreted as complete uncertainty. Generally, the uncertainty is larger in systems with more 
extreme parameters. The dotted red line in each plot shows / se ismic = 10 Hz, the lower bound of sensitivity in Advanced LIGO. The 
dashed green line in each plot shows /2o<y c — the twentieth percentile of the power of the match, defined by Eq. (12). Comparing any plot 
to an accuracy requirement, which may depend on both M toi and <x>h y b, we can extract the maximum sufficient hybridization frequency, 
which suggests the optimal initial orbital frequency. The method for doing this is described in Sec. III. These plots are discussed in more 
detail in Sec. Ill C. 
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For simplicity, the lines shown in the plots use EOB data 
for both waveforms: h a = hf, = /ieob- The line is roughly 
a lower bound for mismatches of MM = 0.20 — the white 
region of each plot. The higher the white region is above 
the /20 % line, the more accurate the PN approximations are. 
At very high frequencies, the white region must approach 
this line simply because the PN approximations break down 
more quickly there. 

D. Mass-dependent target accuracy 

In each of the examples above, the accuracy requirement 
calls for a specific mismatch, regardless of the total mass 
of the system. This is unrealistic for three reasons. First, 
and most simply, the result depends sensitively on a choice 
of mass — in the example of Sec. Ill B, that choice is the 
lower bound of the mass range used in Fig. 4. If we were 
to increase that lower bound to M tot > 50 M Q , rather than 
calling for simulations 87 times longer, we would accept 
the longest current simulations. This sensitivity to the mass 
range shows that it must be considered more carefully than 
an arbitrary choice of plotting range. 

Second, the SNR of an astrophysical signal will depend 
on the total mass. The precise dependence is compli- 
cated, as it involves the shape of the noise curve in Fig. 2. 
However, for masses at which the merger frequency is 
much higher than the detector's low-frequency sensitivity (a 
good approximation for the mass range discussed here), the 
stationary-phase approximation shows that the SNR should 
scale roughly as M^ 6 [39]. If we expect the typical low- 
mass system in real data to have a lower SNR than the 
typical high-mass system, there is no reason to model the 
two with the same precision. More precisely, for optimal 
parameter estimation, the error of a model waveform should 
scale inversely as the square of the SNR of the expected 
signal [32, 40]. 

Finally, the merger rates of real binaries will likely de- 
pend on the total mass simply because formation mecha- 
nisms for such binaries should depend on the total mass [11, 
41^-3]. Though current understanding of such mass depen- 
dence is not great [44], it is an area of active research, and 
will no doubt be improved in the future. In that case, we may 
wish to fold the expected event rate into our target accuracy, 
so that time is not wasted calculating precise waveforms for 
systems we are unlikely to observe. 

For these reasons, useful and efficient accuracy require- 
ments should depend explicitly on M tot . Incorporating 
mass dependence involves reinterpreting the target mis- 



match MM target as a function of mass, rather than as a con- 
stant. For example, we might hope to model a particular sys- 
tem well enough to ensure detection of any binary expected 
in the data, and to allow accurate parameter estimation for 
90 % of binaries. A very crude function that implements 
this idea is 



MM target (M tot ) 



1 

10- 2 

10" 4 ( 



-5/3 



20 M B 



12 < 



M M < 7. 
Mo - 

7 < Mm < jo- 

Mot 

(13) 

The first two cases are inspired by population-synthesis 
results [41] suggesting that basically all equal -mass black- 
hole binaries should have 8 M tot > 7 M Q , and that roughly 
90 % should have M tot > 12 M Q . This function ignores 
binaries with M tot < 7 M Q , and models binaries with 
M tot < 12 M Q just well enough to ensure detection. For 
higher masses, the function scales inversely as the SNR, 
and is normalized to optimally estimate the parameters of a 
20 M Q system having SNR 70. Clearly, more sophisticated 
treatments could incorporate the objectives of detection and 
parameter estimation more smoothly, but this will serve to 
illustrate the idea. 

Regardless of the particular form of the target mismatch, 
we use it to deduce the sufficient value of a»h y b (and hence 
Qq) by plotting the ratio 



max aM {MM (h a ,h b )} 

MM targe t(Mtot) 



(14) 



Where this ratio exceeds 1, the hybridization frequency is 
insufficient. With the target mismatch of Eq. (13), this 
ratio is plotted in Fig. 5. The red curve denotes values 
for which the ratio equals 1. The optimal a> n yb is given 
by the lowest point this curve reaches, GM to tWhyb/c 3 ~ 
0.01 1. Though this frequency is only 47 % higher than the 
frequency deduced in the previous section without incorpo- 
rating mass dependence, it corresponds to a simulation that 
is nearly 3 times shorter — a significant improvement from 
the perspective of the numerical relativist, and hopefully a 
more accurate representation of the accuracy truly required 
for gravitational-wave detector data analysis. This holds 
great significance for systems with large mass ratios, where 
astrophysical considerations suggest that the lower bound 
of M to t w iH b e l ar g e > as mentioned in Sec. IIIC. 



8 Note that several of the papers cited here discuss individual black-hole 
mass or the chirp mass of a binary, rather than the total mass. 

9 A mismatch of 1CT 2 or less ensures detection 97% of the time for a 
signal with SNR at least 7 [40]. However, a real search of detector 
data will use a template bank which does not need to be as accurate 
as this, in the sense that an inaccurate template with different physical 
parameters may happen to match the exact waveform. That is, the error 
of an "effectual" template bank may be far greater than 1CT 2 . The issue 
of detection is complicated, and is discussed at length in Ref . [9] . 
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FIG. 5. Ratio of maximum mismatch to target mismatch. This 
plot shows the ratio between the maximum mismatch between 
various hybrids and a target mismatch given by Eq. (13). For 
values of this ratio greater than 1, the hybridization frequency is 
too high to achieve the target accuracy. The optimal sufficient 
value of tL>h y b is given by the lowest frequency at which the ratio 
is 1 — roughly GM M Whyb/c 3 = 0.011 here. 

E. Mass- and frequency-dependent target accuracy 

Up to this point, we have assumed that the mismatch 
between our hybrids completely describes the uncertainty 
in the final result, after a numerical simulation has been 
done. However, for the stringent accuracy requirements 
quoted above, the numerical simulation will be very long, 
making it difficult to achieve high accuracy in the NR 
portion of the waveform. We may wish to leave a portion of 
the error budget for the NR data, and the rest for the PN data 
and hybrid. Given some understanding of how the NR error 
depends on the length of the simulation, we can incorporate 
that error into our determination of the optimal value of Qq. 
The error of a simulation depends, no doubt, on its length, 
but also on the particular implementation used for that 
simulation and even the computational resources available, 
and is thus beyond the scope of this paper. However, the 
basic idea is a simple extension of the technique discussed 
in the previous section: generalize the target mismatch to 
be a function of both M tot and Wh y b, and plot the ratio 

max aM {MM{h a ,h b )} 

MM targ et(M tot ,W hyb ) 

Again, when this is greater than 1 , a>h y b is too large. For ex- 
ample, the target mismatch may be constructed (crudely) by 
setting an astrophysically motivated target MM Ustro for the 



final waveform, then subtracting the estimated uncertainty 
due to the NR data. Then the permissible mismatch in the 
PN data could be defined as 

MM tar get(M tot , co hyb ) 

:= MM Ustro (M tot ) - MM NR (M tot , w hyb ) . (16) 

This target is crude because the mismatch is not additive, 
but it is a conservative estimate. 

This extension of the method to include dependence on 
o»hyb raises an unfortunate — though realistic — possibility: 
it could be that the ratio of Eq. (15) will never be less than 
1 for masses of interest. For example, if MMNR(M to t, ^h y b) 
increases too quickly as a»h y b decreases, the quantity in 
Eq. (16) will be too small, and thus the ratio in expres- 
sion (15) will be too large. This would indicate that 
the modeling methods, both analytical and numerical, are 
simply too crude to compute the waveform with the desired 
accuracy. 

IV. DISCUSSION 

Figure 4 presents the main results of this paper, showing 
the largest mismatch between any pair of plausible hybrids 
as a function of the frequency of hybridization and the 
total mass of the system. The hybrids are formed using 
the EOB waveform to substitute for the NR waveform. As 
is argued in Sec. I A and shown explicitly for the equal- 
mass nonspinning system in Appendix A, the particular 
choice of substitute does not affect the final results in any 
significant way. The hybrids' inspiral data are supplied by 
TaylorTl-T4 and EOB waveforms, which are attached at 
ajfryb. These hybrids are then scaled to various total masses, 
and mismatches between each pair are calculated. The 
maximum such mismatch is the estimated uncertainty in 
the models. The plots in Fig. 4 assume that the error in 
the numerical portion of the hybrid is negligible, though 
they can be expanded to account for estimated numerical 
errors, as in Sec. HIE. This uncertainty is a reasonable 
proxy for the error — the difference between the model and 
the exact waveform. Given a target uncertainty for the 
complete model, we can deduce the minimum initial orbital 
frequency necessary to achieve that target with a simulation, 
by noting that the relevant value of <x>h y b must be present in 
the simulation data. 

The results show several interesting features. First, the 
uncertainty generally increases as the modeled system be- 
comes more extreme; for a given value of Wh yb an d M tot , 
increasing either the mass ratio or the spin parameter in- 
creases the uncertainty. This is not surprising, since the 
post-Newtonian order of known spin terms is lower than 
the order for non-spin terms [45]. Similarly, PN methods 
are expected to break down for larger mass ratios [37], for 
which more specific methods are necessary [46]. 

More quantitatively, we can relate these results to basic 
accuracy standards for gravitational-wave detectors. To cal- 
ibrate waveforms for detection, accuracies of MM < 0.01 
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are generally called for [32, 40]. Meanwhile, the longest 
current numerical simulations start with G M tot o»h y b/c 3 > 
0.035. The upper-left panel of Fig. 4 shows that hybrids 
created using such simulations would only be sufficient for 
M tot ^ 26 M e . Of course, real detector data is searched for 
a range of system parameters; a real 10 M Q + 10 M Q system 
might be detected by an inaccurate 6 M Q +18 M Q template, 
for example [10]. 10 This dramatically reduces the accuracy 
requirements on a template bank for detection [9]. However, 
while template banks may be subject to loose accuracy 
requirements, numerical relativity will generally be used for 
other purposes — most likely calibrating template banks to 
NR waveforms or hybrids incorporating them. The results 
presented above show that any such calibration is bound to 
exhibit very large errors for low-mass systems unless the 
numerical simulation is very long. 

More stringent demands are placed on waveforms for pa- 
rameter estimation, depending on the SNR of the observed 
signal. For example, modeling the unequal-mass high-spin 
system to high accuracy would require simulating nearly the 
entire in-band signal; the simulation would need to begin 
roughly 40 % above the seismic wall to achieve mismatches 
of MM w 10 -4 . These grim results present discouraging 
prospects for accurate modeling of precessing systems, PN 
approximations for which are known to still-lower order. 

On the other hand, for small values of the mismatch, the 
appropriate value of o»h y b varies almost linearly with M tot . 
The initial orbital frequency Qo required for a simulation 
will then be nearly proportional to the total mass of the 
modeled system, so the length of the simulation will vary 

—8/3 

roughly as M tot . This strong dependence shows clearly 
that target accuracies for a simulation to be used across 
a range of masses should include carefully considered 
mass dependence. Such dependence is incorporated into 
the technique for determining Qo in Sec. HID. For the 
q — 10 systems in Fig. 4, this improves the situation 
dramatically. The smallest black hole observed to date 
has M > 3 M [38], so a binary with q = 10 would 
have M tot ^ 33 M Q , substantially raising the value of Whyb 
required to achieve a given mismatch for astrophysically 
likely sources. In this sense, systems with large mass ratios 
are actually easier to model than comparable -mass systems. 

There are several possible flaws in these uncertainty 
estimates. Most basically, we simply assume that the 
uncertainty in our range of waveforms is a suitable proxy 
for the error in the waveforms. Of course, these models — 
both PN and EOB — may simply be wrong. For example, 
we can imagine that some fundamental error exists in our 



In much of the literature, this difference is highlighted by distinguishing 
between "effectualness" (the match between a given signal and the best 
fit in a template bank) and "faithfulness" (the match between a given 
signal and the particular signal in the template bank with the same 
physical parameters). 



understanding of approximations to Einstein's equations for 
black-hole binaries. In that case, our models may be per- 
fectly precise but entirely inaccurate; the exact waveform 
would lie outside the bounds of our uncertainty estimate. 
Moreover, these estimates depend on the assumption that 
the range of plausible hybrid waveforms is neither too 
narrow nor too broad. The choices made above were based 
largely on coincidences of history; some other reasonable, 
equally accurate, but not-yet-imagined approximant may ex- 
ist, lying far from the approximants used here. Conversely, 
there may be some subtle error in one or more members of 
this group of approximants, leading to unnecessarily large 
uncertainty. Unfortunately, the only obvious way to detect 
such errors is to test the results using very long and accurate 
numerical simulations. 

Taken together, these results indicate that more work 
will be needed to produce accurate waveforms for stellar- 
mass black-hole binaries, even for aligned-spin systems. 
Improvements may come in the form of higher-accuracy 
PN or EOB waveforms, longer numerical simulations, or 
both. This paper has not treated precessing systems simply 
because the production of full waveforms for such systems 
is still in its infancy. No doubt, however, the uncertainties 
are greater than in the cases discussed above. While both 
analytical and numerical relativity have clearly made great 
progress in the past decade, much remains to be done. 
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Appendix A: The (un)importance of the choice of 
ersatz NR waveform 

Section I A presents arguments that the mismatch be- 
tween two hybrid waveforms with ersatz NR data should 
be almost completely insensitive to the phasing of the 
ersatz waveform above the hybridization frequency, and 
only weakly dependent on the amplitude. The key point is 
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FIG. 6. Using real NR data for merger and ringdown. This 
plot reproduces the top-left plot of Fig. 4, using real NR data 
in place of the EOB ersatz NR waveform. The NR data starts 
at frequency a> ~ 0.035, and is extended to lower frequencies 
by hybridizing with a TaylorT4 waveform. The plots are almost 
identical, indicating that the procedure is not strongly sensitive to 
details of the ersatz NR waveform. 
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FIG. 7. Mismatch between the EOB waveform and the NR 
hybrid. This plot shows the mismatch as a function of total mass 
between the EOB waveform used in the body of this paper and the 
NR hybrid used in Fig. 6. The similarity between Figs. 4 and 6 
despite the significant mismatches shown here lead us to conclude 
that the uncertainties shown in those figures are indeed robust with 
respect to the choice of ersatz NR waveform. 



that the two hybrids are identical for oj > a»h y b- It was thus 
argued that the particular choice of ersatz NR waveform 
should not strongly affect the results, as long as the power 
is fairly correctly distributed in the frequency domain. 

One simple way to check this claim is to use a different 
ersatz NR waveform. Here, we will reproduce the crucial 
result of Fig. 4 in the equal-mass nonspinning case with 
different ersatz NR data. In this case, we will substitute 
the EOB waveform with a numerical waveform hybridized 
with TaylorT4 to extend to lower frequencies. The nu- 
merical waveform is the same one introduced in Ref. [47], 
except that Regge-Wheeler-Zerilli wave extraction is used 
to produce h. The waveform is hybridized exactly as in 
Ref. [10]. This hybrid is then substituted for the EOB 
waveform wherever it is called for. 

The results are shown in Fig. 6. Comparing with 
the upper-left plot of Fig. 4, we see excellent agreement 
throughout the plot. The plot shown here does exhibit some 
jagged lines in the range 0.04 < G M tot Whyb/c 3 < 0.045. 
These are evidently due to noise in the waveform itself, 
which appears to be related to junk radiation. That noise 
can easily lead to imperfect hybrids, especially using the 
frequency-alignment scheme of Eq. (8). 

At the very least, this demonstrates that the simplistic 
ringdown-alignment technique used for the EOB waveform 
in this paper (see Sec. B 2) does not significantly affect the 



final results. On the other hand, we might worry that the 
NR hybrid used here is practically identical to the EOB 
waveform used in the main text of this paper, because 
the EOB waveform aligns quite accurately to the very late 
stages of the NR data. In fact, that alignment is misleading, 
because it requires coherence over the relatively short span 
of the numerical data. Judged in terms of the mismatch, the 
NR hybrid and the EOB waveform are quite distinct, shown 
in Fig. 7 as a function of the total mass of the system. 



Appendix B: Details of the implementation 

The results of this paper depend sensitively on accurate 
numerical implementation of the technique. The various 
approximants, their hybrids, and the mismatch must all 
be calculated to high accuracy to ensure that the plots of 
Fig. 4 depict uncertainty in the models, rather than errors 
in the numerical methods. This section outlines the steps 
necessary to obtain accurate results. In short, every attempt 
was made to ensure that the model waveforms were as 
accurate as possible, and each number quoted in this section 
was tested to ensure that making it more stringent had no 
significant effect on the final results. 
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1. Post-Newtonian ingredients 

The TaylorTn approximants used in this paper are based 
on the results of Ref. [26]. In the appendix of that reference, 
the most current and complete PN results are collected and 
expressed in consistent notation. In particular, the orbital 
energy, tidal heating, and gravitational-wave flux are given 
to 3.5-PN order in nonspinning terms, and incomplete 2.5- 
PN order in spinning terms. The spins are assumed to be 
aligned or anti-aligned with the orbital angular momentum. 
These are the basic ingredients to construct the phasing 
of TaylorTl-T4 approximants, as described succinctly in 
Ref. [25]. The orbital phase thus derived is then used in the 
waveform amplitudes described by Ref. [26], which include 
nonspinning terms up to 3-PN order, and spinning terms 
through 2-PN order. 

The main practical concern in constructing these wave- 
forms is producing data on a sufficiently fine grid that accu- 
rate derivatives are available for the alignment step, Eq. (8) 
of the hybridization procedure. For the TaylorTl and T4 
models, this is accomplished by setting a tight tolerance on 
the numerical integration scheme, as discussed in Sec. B 3. 
For TaylorT2, the waveform is evaluated on at least 50 000 
uniformly distributed values of the velocity parameter v, 
which is the independent variable of this model, ranging up 
to v = 1 . Similarly, the TaylorT3 waveform was evaluated 
at 50000 values of the independent time parameter r, 
distributed at roughly uniform intervals of v. Lower values 
of these numbers led to poor hybrids, characterized by noise 
at low values of <x»h y b in the plots of Fig. 4. 

2. EOB model 

The EOB model used for this paper was designed to in- 
corporate recent improvements to the inspiral portion of the 
model, including spin terms, while also remaining robust, 
allowing its application to the somewhat extreme case of 
q = 10, x = 0.95. The primary compromises made in the 
interests of robustness were abandonment of the factorized 
multipolar waveforms of Ref. [48] and coherent attachment 
of the ringdown portion of the waveform. The former 
compromise requires the use of the Pade-expanded flux to 
calculate the phasing of the system, and the standard PN 
multipolar waveforms. These are both reasonable substitu- 
tions: the flux term is used in the EOB code for the LIGO 
Algorithm Library; the PN multipolar waveform should 
still be accurate for most of the inspiral [48]. The latter 
compromise primarily affects the phase of the waveform 
during its very last stages. As was argued in Sec. I A, this 
is unlikely to have any significant effect. In any case, the 
uncertainty of the plausible waveforms is dominated by the 
TaylorTrc approximants in all cases shown in this paper, and 
essentially identical results are obtained when using real 
numerical data for the merger and ringdown, suggesting that 
any error in the EOB model does not affect the results. 

The EOB Hamiltonian used here is roughly the same 
as the one given by Ref. [19], except that nonspinning 



terms in the metric functions A(r) and D(r) are extended 
with new terms from Ref. [33]. Thus, in the nonspinning 
case, the Hamiltonian of this EOB model reduces exactly 
to the Hamiltonian of Ref. [33]; in the spinning case, it 
reduces nearly to the Hamiltonian of Ref. [19]. The angular 
momentum flux is described by Eq. (65) of Ref. [49], where 
the term is given by the Pade expansion of the flux 
from Ref. [26]. The standard formula for u po i e gives very 
poor results for high spins. For this paper, the following 
extension of v po \ e to the spinning case is used: 

6 + 2v 

^Dole = — . - (Bl) 

V(3 + v)(36 - 35v)-at,(8-4v) 

Initial data is set according to Eqs. (4.6) and (4.13) of 
Ref. [17]. Eccentricity is then iteratively reduced to e < 
10 -14 using Eqs. (71) and (73) of Ref. [50]. For the high- 
spin cases, this method does not work directly. Instead, the 
spin is increased in stages. The non-eccentric initial data for 
the given mass ratio is first obtained with x — 0, then used 
as initial data for eccentricity reduction with^ = 0.1. This 
is repeated, incrementing the value of x, until the desired 
spin parameter is reached. That non-eccentric initial data is 
then used to evolve the full inspiral. Reducing eccentricity 
is not only more faithful to the scenario modeled by the 
other approximants, but also allows larger time steps to 
be taken by the numerical integration scheme; significant 
eccentricity would require at least a few steps to be taken 
per orbit. 

The integration ends when the EOB radial parameter is 
smaller than 1, or the radial momentum becomes positive. 
In all cases explored for this paper, the amplitude of the 
resulting waveform reaches a peak, roughly where merger 
is expected, and roughly similar in amplitude to the peak ex- 
pected from numerical simulations. Previously published 
EOB models align a sum of decaying quasinormal modes 
to the rising side of this peak [17, 33, 50]. Those techniques 
do not seem to be sufficiently robust to apply naively to 
the extreme cases discussed in this paper. Moreover, such 
techniques seem to be unnecessary; as argued previously, 
the particular details of the end of the waveform will not 
strongly affect the final results, especially for the small 
portion of the waveform represented by ringdown. For 
these reasons, a simple — though undoubtedly inaccurate — 
method is used to attach a single quasinormal mode to the 
inspiral waveform. The descending side of the amplitude 
peak is used, and the quasinormal mode with the longest 
decay time is attached at the unique point such that the 
amplitude and its first derivative are continuous. 

3. Numerical integration of ODEs 

The TaylorTl, TaylorT4, and EOB waveforms are inte- 
grated numerically by the eighth-order Dormand-Prince 
method implemented in Numerical Recipes [51]. In all 
cases, the absolute tolerance was set to atol = because 
of the vastly different scales of the dependent variables. 
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The value of the relative tolerance rtol was chosen by 
looking at the convergence of the phase of each approx- 
imant. Tolerances from 10~ 4 to 10" 11 gave the same 
results to within small fractions of a radian over the entire 
~ lOOOOOradinspiral. For EOB, rtol = 10~ 6 was chosen 
to be conservative, while still allowing the code to run 
very quickly (less than one second per waveform). For 
TaylorTl and T4, rtol = 10~ 10 was used as a crude 
but effective way of ensuring that output was frequent 
enough to produce smooth derivatives for the alignment 
procedure, Eq. (8). Additionally, dense output was used 
to save 50 intermediate points per time step, which further 
improved the alignment procedure. Practically identical 
results were obtained with the Bulirsch-Stoer integration 
scheme, except that this method could not reliably continue 
into the delicate final few radians of the EOB integration. 
Integration continues until the dependent variables or their 
derivatives reach some unphysical value: for TaylorTl and 
T4, angular frequency is required to remain positive; for 
EOB, radial momentum is required to remain negative, and 
radius greater than 1 . 

4. Fourier transforms and mismatches 

Fourier transforms find two applications in the calcu- 
lation of the mismatch. First, and most obviously, time- 
domain waveforms must be converted to the frequency 
domain for use in the inner product, Eq. (1). Second, the 
match itself is then evaluated by taking an inverse Fourier 
transform. Assuming that the waveforms h a and h/, are 
normalized so that (h a | h a ) = (h/, | h/,) = 1, and combining 
expression (2) with Eq. (3), we see that 



{h a | h b ) = max (h a \ h b ) 

A7\AO 



(B2a) 
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Note that the integral in the last expression is simply the 
inverse Fourier transform of h a h* b /S n (\f\). The maximiza- 
tion over AT involves selecting the largest element (in 
absolute value) of the discrete set produced by the fast 
Fourier transform of that quantity. 

Two concerns drive the application of these Fourier trans- 
forms: aliasing at high frequencies, and Gibbs artifacts at 
low frequencies. To avoid aliasing, the sampling interval 
of the time-domain waveforms must be set by the highest 
frequency of the lowest-mass system of interest. In our case, 
that system has a mass of M tot = 5 M Q . An acceptable 
sampling frequency is half the Advanced LIGO sampling 
frequency: / s = 8192 Hz. On the other hand, avoiding 
Gibbs artifacts at low frequencies requires the waveforms 
to start early enough that the waveform "turns on" outside 
of the LIGO band, and its amplitude is very small at that 
point. Tests with the waveforms used in this paper show 
that an initial frequency of 8 Hz is sufficient to ensure 
accuracy of the mismatch to MM < 10~ 7 . For M tot = 
5 M Q , this corresponds to a dimensionless initial orbital 
frequency of G M tot Q,q/c 3 ~ 1.97 X 10~ 5 . The waveforms 
used in this paper were calculated in dimensionless units, 
starting with that frequency, hybridized as necessary, scaled 
to the appropriate total mass, projected to the positive z 
axis, and interpolated to a uniform time grid with spacing 
At - 1 // s . In extreme cases, these waveforms can consume 
hundreds of megabytes each. Given that five such wave- 
forms need to be compared, and that comparison requires 
significant additional memory, the full memory usage can 
easily reach several gigabytes. Because of the large memory 
requirements, the calculations for Fig. 4 were performed 
on a cluster having ample memory in a single node. To 
use CPU resources efficiently OpenMP [52] was employed, 
which allows very simple alterations of source code to 
incorporate multiple processes — just three additional lines 
of code enabled multiprocessing which resulted in a speed 
improvement by a factor of four. 
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